library(spatialprobit)
dyn.load('test/SpatialProbit.so')
data(Katrina)
#Katrina<-Katrina[1:300,]
Y<-Katrina[,12]
Y[Y==0]<- -1
nb <- knn2nb(knearneigh(cbind(Katrina$lat, Katrina$long), k=20))
listw <- nb2listw(nb, style="W")
W2 <- as(as_dgRMatrix_listw(listw), "matrix")
#W2 <- as(as_dgRMatrix_listw(listw), "CsparseMatrix")
X<-as.matrix(Katrina[,4:11])
X<-cbind(1,X)
W2<-(W2-diag(diag(W2)))

#sarprobit(Y ~ X, W=W2,ndraw=600, burn.in = 100, showProgress=TRUE)




n<-length(Y)
p<-ncol(X)

M<-1000
rho<-matrix(0,nrow=M,ncol=1)
B<-matrix(0,nrow=M,ncol=ncol(X))
res2<-.C("PMCMCspatialProbit",as.double(W2), as.integer(M),as.double(X),as.double(Y),as.integer(n),as.integer(p),beta=as.double(B),rho=as.double(rho))
B2<-matrix(res2$beta,ncol=p)
	
cat("Gibbs")	
M<-50
B<-matrix(0,nrow=M,ncol=ncol(X))
rho<-matrix(0,nrow=M,ncol=1)
res<-.C("GibbsspatialProbit",as.double(W2), as.integer(M),as.double(X),as.double(Y),as.integer(n),as.integer(p),beta=as.double(B),rho=as.double(rho))
B1<-matrix(res$beta,ncol=p)

